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We consider a two-dimensional layer of dipolar particles in the regime of strong dipole moments. 
Here we can describe the system nsing classical methods and determine the crystal structure that 
minimizes the total energy. The dipoles are assumed to be aligned by an external field and we 
consider different orientations of the dipolar moments with respect to the two-dimensional plane of 
motion. We observe that when the orientation angle changes away from perpendicular and towards 
the plane, the crystal structure will change from a hexagonal form to one that has the dipoles sitting 
in equidistant rows, i.e. a striped conhguration. In addition to calculating the crystal unit cell, we 
also consider the phonon spectrum and the speed of sound. As the orientation changes away from 
perpendicular the phonon spectrum develops local minima that are a result of the deformation to 
the crystal structure. 


I. INTRODUCTION 

Quantum degenerate gases of dipolar particles is a fore¬ 
front research topic in cold atomic gases [THIS] • While 
losses can be a problem when dipolar forces are strong 
[3 m m ES], it has been experimentally demonstrated 
that losses can be suppressed by confinement in lower¬ 
dimensional geometries [5]. This opens up the possibil¬ 
ity of realizing interesting lower-dimensional phenomena 
driven by long-range interactions in both few-body [m- 
[23] and many-body [25H55] physics. An overall goal is to 
realize a quantum simulator with long-range interactions 

[asm. 

A recent drive in experiments has been to produce 
dipolar gases using heteronuclear molecules with large 
dipoles moments [Hill ESI- It is then expected that 
one can see strong dipolar effects already before reach¬ 
ing degeneracy. When the dipolar energy scale exceeds 
the energy scale set by the temperature of the system one 
may expect to see crystal formation in the classical sense. 
The regime where strong dipolar forces are the dominant 
feature of the system is the focus of the present paper, 
and we will thus ignore thermal effects. We will be con¬ 
sidering a geometry where the dipolar particles can move 
in two dimensions and where an external field is used to 
align the dipole moments at any fixed angle with respect 
to the plane of motion. We note that near the magic 
angle (to be defined below) where two dipoles on a line 
will have vanishing dipolar interaction energy, a normal 
mode in the crystal will approach zero excitation energy. 
In this regime, we do except quantum and thermal fluctu¬ 
ations to play an important role, but leave this question 
for future studies. 

In constrast to the famous Wigner crystals that have 
been predicted for systems with Coulomb interactions at 
low density [HI ES] , the crystal phases should be found 
at high density when the particles have dipolar interac¬ 
tions. Some earlier works have considered dipoles ori¬ 
ented perpendicular to the motional plane both classi¬ 
cally [46H50] and quantum mechanically [5TH56] . This 
leads to a dipole-dipole force that only depends on the 
relative distance thus strongly simplifying the problem. 


Another line of investigation has been dipoles on a two- 
dimensional lattice [57H5D] . Some of the latter works 
have used different orientations of the dipoles with re¬ 
spect to the lattice. In this case the dipole-dipole in¬ 
teractions no longer cylindrically symmetry in the plane 
but can be highly anisotropic. This may lead to striped 
systems as indicated in different response function ap¬ 
proaches [Ml lfiTH63] . It may as well lead to new unex¬ 
pected properties like the roton minimum in helium [64] . 
We note that similar physical issue can be studied using 
externally oriented magnetic colloids [Ml |66] • 


Naturally, these possibilities complicate both classical 
and quantum mechanical approaches. Classical prop¬ 
erties of crystals are crucial ingredients for subsequent 
quantization EH EH, and investigations of quantum ef¬ 
fects like melting point and heat capacity [HSHTT] . In 
the current work we assume that the dipole moments 
are strong such that classical methods are accurate and 
consider dipolar particles in a plane with no predefined 
lattice, i.e. the particles can move continuously. The 
question then becomes whether, or perhaps more appro¬ 
priately when, the crystal structure starts to change sig¬ 
nificantly from the hexagonal crystal structure that is 
found when the dipoles are aligned perpendicular to the 
plane of motion. 


The paper is organized as follows. In Sec. [IT] we de¬ 
scribe the theoretical model and the parametrization of 
the most general crystal structure in the plane for ar¬ 
bitrary orientation of the dipoles. Sec. jlllj outlines the 
results for the crystal structure itself. In Sec. jlVj we 
elaborate on the properties of the minimal energy crystal 
configuration. We calculate the phonon spectrum by out¬ 
lining a proper parametrization of the reciprocal lattice 
before presenting the spectra along with the sound ve¬ 
locity. Sec. [V] contains a short summary and an outlook. 
The technical details of the calculations of the phonon 
spectrum are presented in Appendix [Aj 
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FIG. 1: The angular coordinates, 9 and 4>, describing the 
direction of the dipole moment, rh. 


FIG. 2: Illustration of two identical aligned dipoles on the x- 
axis separated a distance r from each other. Both directions 
are described by the same {9, (f>) as in Fig. 


II. CRYSTAL PARAMETRIZATIONS 

In order to determine the crystal structures for dif¬ 
ferent aligned dipole moments of the particles, we need 
to develop a convenient parametrization of the geome¬ 
try. We imagine a system with a number of particles 
distributed on a section of a plane in lattice-like struc¬ 
tures. Instead of the number of particles in the allowed 
space it is much better to use the density of particles 
(an inverse area in 2D) as a decisive parameter. The 
anticipated regular structure can then easily be broken 
down into unit cells repeating each other across the plane. 
The choice of parametrization is intrinsically connected 
to the interactions between the particles. We therefore 
first examine the details of the dipole-dipole potential 
and subsequently in the next subsection we describe the 
parametrization of the unit cell. 


A. Dipole-dipole interaction 

All dipoles are assumed to be identical with a given 
dipole moment of size D. They all have the same orien¬ 
tation induced by the application of an external magnetic 
or electric field. The direction of the dipole moments is 
given by the angles (0, (j)) as illustrated in Fig. The 
dipole moment itself, rh is then described by 

/sin(0) cos{4>)\ 

rh = D I sin(0) sin(())) ] . (1) 

V cos(6») J 

The interaction energy between two of these dipoles, rhi 
and m 2 , is 


I = C 


f rhi ■ m 2 3 (mi • r) (m 2 • r) 

\ fp3 rp5 


( 2 ) 


where r is the vector between the two dipoles, r = |r], and 
C = ^ and C = for magnetic and electric dipoles. 


respectively. Here /Tq and eg are the vacuum permeability 
and vacuum permittivity, respectively. Below we will as¬ 
sume magnetic dipoles. Results for electric dipoles can be 
obtained by simply making the substitution /rg —>■ l/eg. 

We now assume that the two dipoles are placed on the 
cc-axis at a distance r and described by the angles 0 and 
4>. For identical dipole moments (mi = m 2 ) = rh the 
potential energy in Fig. is reduced to 


(3) 

since r is taken along the a;-axis. We thus see that the 
interaction energy is positive as long as 

1 — 3sin^ 0cos^ 0 > 0. (4) 

We therefore conclude that the interaction is repulsive 
for all cj) when 6<9c = sin“^ (l/-\/3) = 0.61548. Thus, 
collapse due to attractive inverse cubic interactions can 
not occur for 9 < 0.61548. We note that the sin^ 9 depen¬ 
dence implies a symmetry around 9 = t: 12 corresponding 
to reflection in the xy-plane where all dipoles are assumed 
to be located. 


CD^ 


(1 — 3sin^ 0cos^ (j)) , 


7 = c(?-hh, 


B. Unit cell 

When comparing different crystal structures we must 
consider the conditions we impose on the system. We 
imagine that the crystal forms from a gas phase that con¬ 
denses by cooling. Consequently the crystal structures 
most likely prefer to have a constant density of dipoles 
in the plane, i.e. we do not have clumping or holes in the 
crystal. We thus assume an ideal uniform crystal struc¬ 
ture. It is of course very interesting to study the effect of 
crystal defects in both a classical and quantum mechan¬ 
ical setting but this is beyond the scope of the present 
study. 
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FIG. 3: Illustration of the parametrization used for the unit 
cell in the crystal xy-plane where the 2 -axis is perpendicular 
to the plane with the particles. The unit cell is shown in 
the bottom left corner with bold lines. The parametrization 
contains two variables, 9' and b. 6' determines the angle in the 
lattice, a and ab are the distance between nearest neighbors 
in the x- and the ^'-directions, respectively. With this choice 
we can parametrize all uniform crystal structures. 

We therefore aim at determining the preferred regular 
crystal structure for a given alignment of the many dipole 
moments distributed to produce a given average density 
in the xy-plane. We then need to find the configuration 
that minimizes the total interaction energy for such spec- 
ihcation of the system. The total interaction energy, 17, 
is 

U=\'^^I{R-R'), (5) 

RR' 

R^R' 

where R and R' are position vectors of each pair of dipoles 
interacting through the potential, /(r), in Eq. ([^. 

We need a sensible unit cell that can be used to 
parametrize the crystal structure for given external align¬ 
ment. An important limit to consider is the one of per¬ 
pendicular dipoles, i.e. when all the dipole moments are 
aligned perpendicular (9 = 0) to the plane containing all 
the particles. In this case the dipole-dipole interaction is 
purely repulsive with a r~^ behavior that only depends 
on the relative distance. The system will try to maxi¬ 
mize the distance between each pair of dipoles and not 
allow any clusterization or irregular correlations (at least 
in the ideal crystal we investigate). Finding the min¬ 
imum energy configuration for this system is therefore 
the same problem as the packing of spheres. This prob¬ 
lem is known to have a hexagonal solution in 2D |72)and 
we therefore must be able to capture this geometry with 
our unit cell parametrization. 

To allow sufficient flexibility we choose the unit cell 
shown in Fig. The two variables 9’ and b allow us 
to describe all uniform crystal structures with a as the 
crystal lattice length. Symmetry allows a restriction of 
0’ to the interval [0,7r/2[ and b to be less than unity, i.e. 
]0,1]. The value of a can be determined as function of 0’ 
and b to maintain a given average density of dipoles in 
the resulting crystal. This can be imposed by insisting 
that the unit cell area remains constant as we change 
the parameters 9' and b. The density, a, is inversely 
proportional to the area, a^b sin 9', of the unit cell. The 


constraint on a is then 

The two unit cell basis vectors, at, connecting nearest 
neighbors in the two principal directions are 

^ (a\ ^ (ab cos 9'\ 

ai=(^oj’“2=(,a&sin0'j' 

We may now exploit the symmetries of the unit cell. As 
mentioned before there is a symmetry around 9 = 7r/2. 
There is also a symmetry resulting in U {(p) = U {p tt) 
where the energy is invariant under translation of p by 
TT. This can easily be seen by the invariance of rotating 
the unit cell 180 degrees around the z-axis, e.g. rotating 
around the central point in Fig. It can also be seen 
by looking at the expression for the dipole interaction 
in Eq. <§■ The dependence on (p goes as cos^ (p and 
this is invariant under rotations by tt. We can therefore 
reduce the parameter space and work with 9 € [0,7r/2] 
and (p e [0, tt]. 

It is important to note that the external alignment of 
the dipoles fixes 9 as the angle between the z-axis and 
the direction of the external electric or magnetic field. 
The other angle, (p, describes the angle between the ar¬ 
bitrarily chosen cc-axis of the unit cell and the projection 
on the crystal plane of the aligned dipole moments. The 
calculational procedure will therefore be to fix 9 and cp, 
and then minimize the energy with respect to 9' and b. 
This will for each 9 produce an energy as a function of (p 
and the minimum will then define the preferred crystal 
structure for this given 9. 

In general a complicated topology may invalidate this 
procedure of fixing one parameter, minimizing with re¬ 
spect to two other parameters, and afterwards find mini¬ 
mum of the resulting function of the first parameter. All 
three variables should perhaps be varied simultaneously 
to find true global or local minima. To confirm that our 
procedure provided correct minima we tested by numer¬ 
ical variation of the parameters around the minima. 

III. CRYSTAL STRUCTURE 

The structure of uniform crystals with constant pla¬ 
nar density of dipoles has been determined as function 
of polarization angles. We shall here present the result¬ 
ing configurations for a number of angles that illustrate 
the general behavior. We will start with 9 = 0 (indepen¬ 
dent of (p). Subsequently we increase 9 up to the critical 
value, 9c, for collapse due to unhindered small-distance 
attraction. We shall investigate how the crystal evolves 
for different external alignments of the dipoles, that is for 
the interval 0 < 9 < 9c where (p, 9' and b assume corre¬ 
lated values minimizing the total energy. We emphasize 
again that 9' and b describe the crystal structure in Fig.[^ 
whereas (p describe the orientation of the crystal in the 
plane with respect to the external polarization direction. 
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FIG. 4: The total energy for 0 = 0 (independent of <(>) as a 
function of b and 6' is here shown in units of Uo = . 

The constant two-dimensional density, a, is an inverse area 
and consequently is an inverse cubic length. Only a 

section of the allowed values of b and 9' is shown as it is easily 
seen that the energy increases rapidly beyond what is shown. 
The lowest values are dark blue and increasing energies follow 
the colors of the rainbow. 



FIG. 5: The hexagonal lattice structure for 6 = 0 with the 
specified crystal parameters, 9' and b, as defined in Fig. 


A. Perpendicular dipoles: 0 = 0 

First we consider the structures arising for 0 = 0 where 
the dependence on (j) drops out. It is indeed the simplest 
case yet it is a relevant benchmark to demonstrate the 
validity of our method. The total interaction energy from 
Eq. @ as a function of 9' and b is shown in Fig. En¬ 
ergies are plotted in units of Uq = where the 

density of dipoles, cr, is kept constant in all calculations. 

As expected, we find that the minimum in energy oc¬ 
curs when O' = tt/S and 6 = 1. The unit cell of the 
corresponding hexagonal lattice structure of lowest en¬ 
ergy is shown in Fig. Note that for 0 = 0 there is a 
triangular symmetry which is the result of the rotational 
symmetry in the plane for 0 = 0 which is broken down to 


the point-group symmetry of the hexagonal lattice when 
the crystal forms. 

A rather flat valley extends in Fig. from the min¬ 
imum in the direction of smaller 6 for fixed 9' = tt/S. 
The walls increase rather steeply on both sides of this 
valley. The crystal is therefore relatively soft towards 
moving the upper (and lower) dipoles in Fig. further 
away from (or closer to) each other and the central line 
of dipoles along the direction of 9' = tt/S. We emphasize 
that precisely the same unit cell is obtained for 0' = 27r/3. 
Consequently a similar minimum would show up by ex¬ 
tending Fig. |4]to larger values of 0'. The barrier between 
the two minima is already indicated and almost reached. 
No other direction parametrized by 9' provide the same 
structure for the same density. 

Small changes of 0 away from zero would break the 
hexagonal lattice symmetry and introduce a preferred 
direction connected to the direction of the external field 
that aligns the dipoles. The soft directions along 0' = tt/S 
and 0' = 27r/3 then must be either followed or overpow¬ 
ered. The hexagonal structure must be distorted or com¬ 
pletely changed. 


B. Tilted dipoles: 0.1 < 0 < 0.6 

We now move the polarization direction, 0, away from 
being perpendicular to the crystal plane. The crystal has 
to adjust by finding the minimum energy with respect to 
6, 0', and (j)- For convenience we minimize with respect 
to b and 9' for each value of (p. The resulting functions 
are shown in Fig. for three different values of 0. 

The energy for the smallest polarization angle, 0 = 0.1, 
is exhibited in the upper panel. The variation with p is 
exceedingly small with a decrease from the 0 = 0 result 
of 0.7072 by around 1%. 

The case of small angle polarization is nevertheless the 
first step away from the highest symmetry and towards 
collapse of the crystal structure. The case of 0 = 0.1 
in Fig. 1^ already exhibits all the tendencies to develop 
pronounced minima in the potential energy surface. We 
observe the global minimum at (/) = 1.083, which corre¬ 
sponds to 9' = 1.056 and 6 = 0.99. The structure for this 
minimum is then visualized as a polarization direction, 
(/), tilted in the direction of 0' = 7r/3. 

We also note additional (local) minima around 4> = 
2.058 and p = 0.542. The first of these corresponds to a 
structure with p = 9' = ‘2 'k that is the same geometry 
as the global minimum but for a different unit cell. The 
second minimum corresponds to « tt/B, 9' « 7r/3, and 
6 « 0.99, that is an apparently metastable structure. We 
emphasize that all minima are tested to be real minima 
by direct computation of the energy surface in a cubic 
grid around the minimum points. In any case the energy 
differences are very small and sometimes at the limit of 
being significant. 

Still it is of interest to understand the emerging struc¬ 
tures. The minimum energy configuration for p = 0.542 
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FIG. 6: The interaction energy in Eq. ([^ in units of Uo de¬ 
fined in Fig. 1^ for 6 = 0.1,0.15,0.35 (upper, middle, lower 
panels) as a function of 4> after minimization with respect to 
6' and b. 

is the same as for the two deeper minima but now the 
polarization direction is half way between those of the 
preferred global minimum, cj) « O' = tt/S, and the ex¬ 
treme linear structure, (() « 0' = 0. The symmetry is 
very high as the polarization points directly to towards 
the opposite, furthest away, dipole in the unit cell. 

Going away from the flat region for intermediate (j)- 
values, both towards smaller and larger values of (j) we 
find a steeply increasing energy. The curve beyond (/> = tt 


FIG. 7: 9' and ij) is here shown for the global minima as a 
function of 9. 


continues precisely as for ^ = 0. The maximum is rela¬ 
tively high and associated with a structure corresponding 
to the hexagonal structure which is energetically unfa¬ 
vored. 

Tilting the dipole to 0 = 0.15 leads to larger variation 
of the potential energy as shown in the middle part of 
Fig.§ The deepest minimum associated with the op¬ 
timal crystal structure still appears for about the same 
structure values, (j) = 2.058, O' = 1.025 and 6 = 1, as 
for 0 — 0.1. The minimum with the same crystal struc¬ 
ture for about </> = 0.975, 9' = 1.067 and h = 0.98 is 
fortunately also found with the same energy. The over¬ 
all features in the energy as function (j) remain the same 
as for 0 = 0.1 but now determined with better relative 
accuracy. The local minimum at about (j) « 0.542 has 
almost disappeared. 

Proceeding to the larger polarization angle of 9 = 0.35 
we find the energy shown in the lower part of Fig. The 
two optimal minima are now rather pronounced at the 
same (/(-values of ^ = 2.028 and 1.133. The two local 
minima at about (j) = 0.600 and 2.500 are now clearly 
seen in the figure. They correspond to crystal structures 
described by 0' « 1 and b « 0.35. These local minima 
are sensitive to the influence from dipoles at very large 
distances which tend to increase the minimum values and 
perhaps eventually wipe them out completely. 

The minimum energy configurations for the different 
0-values all correspond to 5 « 1. The variations of O' and 
(j) as functions of 9 are shown in Fig. The similarity 
is striking as the two equivalent minima correspond to 
(p K, O' and (j) « 20'. As already mentioned this is not 
surprising as the same structure arises and the polariza¬ 
tion direction points along the highest symmetry axis of 
the crystal structure. 

The two identical optimal structure changes systemat¬ 
ically with 0 as illustrated in Fig. It is seen clearly 
that the structure evolves continuously from the hexag¬ 
onal structure in Fig. and towards lines of well sepa¬ 
rated dipoles. The ultimate linear configurations place 
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FIG. 8: The structure for d = 0,0.2,0.45 and 0.6. The color 
fits with the structures shown. 


the dipoles at pairwise minimum energy positions with 
respect to each other. 

The change in structure seems to be accelerating as 
Q becomes larger and approaching the structure where 
each row of dipoles would collapse on itself after crossing 
the threshold into the range of inverse cubic attraction at 
zero distance. This configuration is very near the point 
Qa = 0.615 where attraction between two dipoles would 
occur, and it would seem logical that this is the point 
where the structure would collapse. The last configura¬ 
tion we investigated was for Q — 0.6 where the angle O' 
is approaching the dive down towards zero as seen in in 
Fig.0 

The overall change in the system makes sense intu¬ 
itively. As the dipole-interaction weakens as 0 grows it 
makes sense that the dipoles would try to place them¬ 
selves behind each other at the minimum energy posi¬ 
tions with the nearest neighbours. The more it weakens 
the more favourable this positioning becomes and the 
structure will change accordingly. 


IV. PHONON SPECTRA 

When the structure has been determined for a config¬ 
uration of the dipoles it will be of interest to determine 
the phonon spectrum. This requires calculations of the 
frequencies of the normal modes of the lattice vibrations 
around stable minimum configurations. We shall first 
indicate the theoretical derivation and provide expres¬ 
sions and calculational procedure for spectra and corre¬ 
sponding speed of sound. The details can be found in 
Appendix Second we exhibit the symmetries of the 
reciprocal lattice as function of O' and b. Finally we dis¬ 
cuss the calculated frequencies, for different polarization 
angles, 0, as functions of wave number in the two normal 
mode directions. 


A. Formulation 

The interaction energy is calculated by use of Eqs. (H) 
and ([^ for a given equilibrium structure where the 
dipoles are located at a series of lattice points, Rq. We 
move the dipoles a small distance away from their re¬ 
spective equilibrium positions, u{Ro), such that R = 
Ro + u{Ro)- Then change of the pairwise distance be¬ 
comes Su = u{Ro) — u{R'q) which assumed small allows 
expansion of the energy in Eq. ([^ to second order in Su. 
The zeroth order term will not contribute to the dynam¬ 
ics. The first order term will vanish as we expand around 
minimum positions. 

We now assume equilibrium at R and omit from now 
on the index, "0". The gradient of the second order term 
provides the force on the amplitude, u, and the equation 
of motion becomes 

M^{R) = = -'£mR-R')u{R'), (8) 

R' 


where M is the mass of one dipole, and the elements of 
the Z?-matrix are defined by 


D^,{R-R') 


d^I 

dUfj,{R)du^{R') 


(9) 


We search for periodic solutions, u{R), in a direction de¬ 
scribed by the unit vector, e, by inserting the ansatz, 

uiR) = (10) 


into Eq. ([^ for a constant given wave vector, k. The 
resulting two-dimensional eigenvalue equation gives us 
two solutions corresponding to vibrational frequencies, 
uj{k), and related to the two normal mode directions, e, 
parallel and perpendicular to the wave vector, k. 

From the eigen frequencies determined for the central 
symmetry point, F in Fig. we calculate the correspond- 
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ing sound velocities, 


where di are the lattice vectors in Eq.Q. This leads to: 


duj{k) 

dk 


k=0 


( 11 ) 


bi = 


277 


^ a tan 9' 


ho = 


0 

277 

> ab sin 9' 


(13) 


for the two normal modes, transverse and longitudinal 
waves. 




FIG. 9: The reciprocal lattices for 9 = 0 corresponding to 
9' = 7r/3 and 6=1 (upper panel), and 9' > 7r/3 and/or 6 > 1 
(lower panel. The lines represent the first Brillouin zone for 
each particle. Symmetry points at F, X and J are observed. 
These are respectively the point between two dipoles and the 
corner between three Brillouin zones. 


B. The Reciprocal Lattice 

To analyze the phonon spectrum we need to know the 
points of interest in the reciprocal lattice. We therefore 
calculate the reciprocal lattice as a function of 9' and b. 
We define 6 i an 62 as the reciprocal lattice vectors. From 
the definition of reciprocal lattice vectors we have 

O'i ■ bj = 6ij2j7T^ (f^) 


The symmetric structure oi 9' = ^ and 6=1 shown in 
Fig. [5] then corresponds to the reciprocal lattice in the 
upper part of Fig. Three different symmetry points 
are immediately recognized and marked by T, X and J 
in the figure. They are respectively the point (0,0), the 
point between the two dipoles at ( 0 , 0 ) and 61 + 62 , and 
the corner between three Brillouin zones for the dipoles 
at ( 0 , 0 ), 61 and 61 +62- 

Increasing the value of 9' above tt/S moves the 61- 
vector closer to the x-axis. This results in a twisting 
of the reciprocal lattice as shown in lower part of Fig. 
The symmetry points, T, X and J, for the perpendicular 
dipoles solution are also present for all other values of 
6 and 9'. Decreasing 9' would lead to a twisting of the 
reciprocal lattice in the opposite direction of going from 
upper to lower part of Fig. 

These symmetry points will be natural target points for 
the wave vector in investigations of the phonon spectra 
of the different structures. 


C. Results for perpendicular dipoles: 9 = 0 

The phonon spectra are calculated along the closed 
path passing through symmetry points T, X and J shown 
in upper panel on Fig. The wave vector, k, always 
begins at the origin, T. First it increases in size in the 
direction of X, then it follows the line towards J, and 
finally it maintains this direction while decreasing size 
until zero at the starting point, T. 

The phonon spectrum is then calculated and shown in 
Fig. for perpendicular dipoles of 0 = 0. Two fre¬ 
quencies arise from diagonalizing the 2 x 2 D-matrix. 
The two types of vibrations are denoted longitudinal and 
transverse corresponding to small amplitude motion of 
all dipoles along and perpendicular to the fc-direction, 
respectively. 

The two computed frequencies both increase from zero 
at the symmetry point, T. The longitudinal mode show 
the fastest increase of energy, reaching a maximum at the 
edge of the Brillouin zone, when the wave vector reaches 
its maximum at the symmetry point X. The transverse 
frequency increases slower but with a similar flat region 
before reaching X. 

With the wave vector end-point continuing from X to¬ 
wards J we find a smoothly decreasing longitudinal fre¬ 
quency and a more abruptly, yet continuously, increasing 
transversal frequency. 

As the symmetry point, J, for the hexagonal lattice is 
approached the two frequencies meet each other in one 
degenerate value. As this is a highly symmetrical point 




























FIG. 10: The phonon frequency in units of loq = \J 
for perpendicular dipoles of 0 = 0 as function of the path 
followed by the wave vector, k. The direction of k corresponds 
to the path with end-points marked in Fig. that is from F 
over X, towards J and back again to F. The relative sizes on 
the x-axis are arbitrary. 

for the hexagonal lattice such a behavior would be ex¬ 
pected here but it is also only expected for this particular 
case. 

The two frequency curves continue smoothly through 
this point, but interchange correspondingly vibrational 
character between transversal and longitudinal mode. 
The longitudinal mode continues to reach a maximum 
before decreasing towards zero at T, while the transver¬ 
sal mode decreases directly towards zero at the central 
point. 

Both frequencies approach zero at the symmetry point, 
r, although with different rates. Thus no optical branch 
is found in agreement with what is expected for a single 
particle basis. 

D. Results for tilted dipoles: 0.1 < 0 < 0.6 

We now slowly increase 9 above zero. For small values 
the structure is rather similar to that of d = 0. Ac¬ 
cordingly the changes in the phonon spectrum is small 
although noticeable. The complete hexagonal symmetry 
is broken and a ^-dependence of the interaction energy 
appear for finite values of 9, see Fig. The symmetry 
is no longer perfect even for a small value of 9 = 0.1, 
but the symmetry points remains in the reciprocal lat¬ 
tice sketched in Fig. The degeneracy in the point J 
is lifted and the two branches begin to move away from 
each other at this point. 

Increasing 9 changes the structure substantially espe¬ 
cially towards the critical value 6*c- The results for in¬ 
creasing 9 are shown in Fig. [m for three larger illus¬ 
trative values. For 9 — 0.25 the degeneracy at point J 
has now clearly disappeared and a significant gap has ap¬ 
peared. The spectrum still maintain the same features 



FIG. 11: The phonon spectra for tilted dipoles with polar¬ 
ization directions corresponding to 9 = 0.25 (upper panel), 
9 = 0.40 (middle panel), and 9 = 0.55 (lower panel). 


and the same overall appearance. However, now a small 
minimum appears on the longitudinal branch close to the 
point J. This minimum is present for all spectra in the 
interval from 9 = 0.2 — 0.35. After appearance it first 
becomes deeper but as 9 grows so does the frequency at 
J and for 9 = 0.35 the minimum is altogether washed 
out. The significance of this minimum could possible be 
a signal of roton dynamics in the crystal in analogy to 
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FIG. 12: 


•J 


M 


: The speed of sound as a function of 9. 

3/2 


Here uo = 


that of helium |64j . Such rotons have been predicted and 
discussed also for particles with dipolar interactions [731- 
UM- It is clearly related to deformation of the crystal, 
but we postpone more detailed investigations to future 
work. 

Increasing 9 to 0.4 and beyond, we see that the spec¬ 
tra in Fig. [m change appearance. The point X moves 
closer to F in Fig. We note that no new minima ap¬ 
pear beside the kink at J in the transversal frequency for 
6 — 0.4. The spectrum for 9 = 0.55 is apparently more 
distorted with kinks or abruptly changing values. How¬ 
ever, we emphasize that the a;-axes on Fig. do not 
reflect real distances. It is merely the result of the choice 
of discrete wave vectors k, and the subsequent distance 
between points on the figures. Therefore the variation 
is not a realistic signal representative for any physical 
effect. 


E. Sound velocities 


With the phonon spectrum available we calculate the 
speed of sound, v, from Eq. 0 in the point, F, in both 
transverse and longitudinal directions. The results are 
seen to be the derivatives of the frequency curves in Figs. 
[T0|and[TT| Each configuration then has longitudinal and 
transversal values for the speed of sound as shown in in 


Eig. 12 as functions of 9. For the hexagonal structure, 
d = 0, we find v = 1.19481uo and v = 0.364272vo, re¬ 
spectively, in the natural units of vq = \J ^ 

Both longitudinal and transversal velocities decrease 
continuously and relatively slowly with 9. There is a ten¬ 
dency to speed up the reduction as the critical polariza¬ 
tion angle is approached. This is intuitively understand¬ 
able, since the structure for 9 ^ 9c approaches separated 
straight lines as seen in Fig. The frequencies decrease 
when the instability is approached, and correspondingly 


the speed of sound must decrease. This instability in also 
clear to see in bottom panel of Fig. [m where we see the 
normal mode going towards zero energy at the point X. 


V. SUMMARY AND OUTLOOK 


Using a classical analysis we have determined the struc¬ 
ture, the phonon spectrum, and the speed of sound for 
dipolar crystals in two dimensions with different polar¬ 
ization of the dipole moments with respect to the two- 
dimensional plane of motion. This analysis has been re¬ 
stricted to polarization angles less than the critical angle 
for dipole-dipole collapse. This ensures that the overall 
interaction is repulsive so that one avoids collapse due to 
head-to-tail attraction of the dipoles. 

For perpendicular orientation, we find the expected 
hexagonal crystal structure. As the polarization angle 
increases, we observe a distortion of the crystal lattice 
and the structure changes towards a more square type of 
lattice. This is associated with the increasing tendency 
for the system to prefer stripes of dipoles. The physics 
of this deformation is associated with the decreasing re¬ 
pulsive energy felt by two dipoles as they are tilted away 
from perpendicular orientation. The optimal energy con¬ 
figuration therefore becomes one where dipoles are placed 
on lines and effectively makes a striped system. Indica¬ 
tions of stripe formation have also been found in different 
response function approaches [6TII63j . 

As we have demonstrated above, the transition from 
hexagonal toward a striped character of the system is par¬ 
ticularly clear in the phonon spectrum and in the speed 
of sound of the system. Measurements of phonon disper¬ 
sions and speeds of sound would therefore be a promising 
way to experimentally probe the crystal. Interestingly, 
we find that the phonon spectrum for increasing tilt an¬ 
gles develops some local minima analogous to the roton 
minima seen in Helium. These new minima occur in po¬ 
sitions different from the ones found in the case of per¬ 
pendicular orientation and are thus associated with the 
deformation of the crystal structure. 

Our study provides the basis for including quantum 
mechanical effects in the system. This becomes impor¬ 
tant once the dipolar interaction strength is compara¬ 
ble to the kinetic energy of the zero-point motion in the 
crystal which occurs for weak dipole moments. When 
the zero-point energy is large we expect strong quantum 
fluctuations and suppression of crystal formation. In¬ 
cluding these effects, one can investigate melting transi¬ 
tions and heat capacity due to quantum motion. This can 
be done by using canonical quantization of the phonon 
modes and/or via a Lindemann criterion approach. This 
will be the subject of future work. 
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Appendix A: Derivation of the Phonon Spectrum 

The derivation will follow the derivation in chapter 22 
in Ashcroft and Mermin (1976) [72]. 

As mentioned before the energy of a structure is cal¬ 
culated as 

i=\'£4>{R-R')- (A1) 

RR' 

R^R' 

Here R and R' are the position vectors of each dipole in 
the structure and (p is the interaction potential between 
two dipoles. We introduce a small perturbation in the 
system as 


where (pf^u^R - R') = ^ ^ (R - R')- Because 

4>fj.i^iR-R') = 4’tJ.i^iR'-R) ^ Ufj,{R)4>fj.„{R-R)uv{R) = 
u^{R')(j)fj,v{R — R')u^{R'), and because R and R' ex¬ 
press the same points it will, per symmetry, be that 
u^{R)4>^v{R - R')uv{R') = u^{R')(j)f^v{R - R')u„{R). 
The harmonic term will again be rewritten as 

^ U^{R)P^^R - R')MR) 

RR' / A o^ 

R^R' 

- Uf_i{R)(j)f^^{R - R')uy{R'). 

A new notation is introduced so only one term is present: 


r{R) = R + u{R). 

Insert this in the expression for the energy. 


(A2) 


^ U^{R)D ^,{R - R')u,{R'), (A9) 

RR' 


^ ^ E (AS) 

RR' 

R^R' 

= b E - “(^'))- (A4) 

RR' 

R=jtR' 

An expansion of the potential is made in u{R) — u{R'), 

f{f+a) = /(fO + a • V/(r) -h ^(a • V)V(fO ''' (A5) 

All terms of order 0{u^) are removed as a harmonic ap¬ 
proximation a is made. This is appropriate as vibrations 
should be small for zero Kelvin. The zeroth order term 
will in principle be an infinite term but it will not be rel¬ 
evant for motion. The linear term will not be relevant as 
we assume all particles to be at minimum positions and 
will therefore be zero. The Harmonic term, 

^ X[U,{R)-U,{&)] 

RR' 

R^R' 


can be rewritten to something more useful 

E {Ut.{R)<Pf.AR - R>AR) 

RR' 

fi^L'—x,y 

R^R' 

+ Uf_,{R')4>fj.^{R - R')u^{R') 

- Uf_i{R)(j)f^^{R - R')uy{R') 

- Uy^{R')Py.u{R - R')u^{R)), 


(A6) 


d^piR - i?0 

drfj,dr„ 


u^{R) - u^{R')]. 


where Dv^{R-R') = 4>vy.{R- R") - PyiuiR- 

R') and the sum allows R = R'. Introduce a matrix 
notation instead: 

jHarm = 1 ^ u{R)Il{R - R')u{R'). (AIO) 

RR' 

This matrix is symmetric as a result of ^ tr^dr^ ^ “ 
a^4>{R-R') 

drydr^ 

It it also clear that = [R — R') can be expressed 
as: 


D^,{R-R') 


d^I 

du^{R)duu{R') 


(All) 


Now two equations of motion can be made for each par¬ 
ticle 


pi pHarm 

Mu^{R) =F^ = ^ = - E ^/^-(^ - R)MR')- 

d^,iR) 

(A12) 

Here the differential has removed two sums in the expres¬ 
sion, M is the mass of each dipole and is the force in 
the ^’th direction. These equations can then be collected 
in vector notation as 

Mu{R) = -J2mR- R)uiR')- (A13) 

R' 


Solving these equations require us to make the ansatz: 


u{R,t) = 


(A14) 
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where k is a wavevector. This results in: 

= - R')u{R',t) 

R' 

= - ^ a{R - 

R' 

= - ee"*“* ^(R - 

R' 

Muj‘^e=eY - R')e"^'^^’~^^ 

R' 

=eYUiR' - 

R' 


(A15) 

(A16) 

(A17) 

(A18) 


where D{k) = ^r" D{R")e = 

—2J2ii»D{R'')sin'^{^k ■ R"). The expression has 
now been reduced to a eigenvalue equation which can be 
solved for a given k. If the structure would be of finite 
proportions there would be restrictions on the allowed 
values of k but if infinite it could have any value in the 
first Brillouin zone. As the matrix has real values, is 
symmetric and of dimension 2 there will be two real 
eigenvalues with two orthogonal eigenvectors which 
fulfils: 


(A19) 

(A20) ^{k)il{k) = \i(k)el(k). (A23) 


(A21) 


From this Wi(fc) belonging to €i{k) can be found as 


Because both R and R' describe all positions in the infi¬ 
nite structure we now make the substitution R” = R — R' 
and get: 

Muj‘^e = eY = ^(fc)e, (A22) 

R" 
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